Pour ce TP, nous allons utiliser les librairies de R suivantes :
library(reticulate)
library(ggplot2)
library(gridExtra)
library(plotly)
library(FactoMineR)
library(factoextra)
library(dplyr)
library(klaR)
library(reshape2)
library(mclust)
library(Rmixmod)
library(cluster)
et les librairies de Python suivantes :
import numpy as np
import pandas as pd
import plotly.express as px
import matplotlib.pyplot as plt
import prince
from kmodes.kmodes import KModes
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
from yellowbrick.cluster import SilhouetteVisualizer
import gower
from scipy.spatial.distance import squareform
from scipy.cluster.hierarchy import dendrogram, linkage
from sklearn.cluster import AgglomerativeClustering
from sklearn.metrics.cluster import adjusted_rand_score
L’objectif de ce TP est d’étudier quelques stratégies pour faire du clustering d’individus décrits par des variables qualitatives.
Nous allons nous intéresser au jeu de données de champignons disponible sur le site de l’UCI : mushroom. Pour limiter le temps de calcul et faciliter l’interprétation durant le TP, nous allons considérer que \(800\) individus (\(400\) de chaque classe) parmi les \(8124\) disponibles et l’on se restreint aux variables suivantes :
Dans la suite, on note la matrice des données \[\mathbf{X}=(x_{ij})_{\begin{array}{l} i=1,\ldots,n\\ j=1,\ldots,p \end{array}}\] où \(x_{ij}\) est la modalité prise par le \(i\)-ème individu pour la \(j\)-ème variable qualitative.
Question : Chargez le jeu de données disponible dans
le fichier mushroomTP.csv déposé sur la page moodle du
cours. Conservez à part la variable Class.
Data <- read.table("...")
Classep <- ...
Data <- .... # jeu de données sans la variable Class
Correction :
Data <- read.table("mushroomTP.csv", header = T, sep = ",", stringsAsFactors = T)
classesp <- Data[, 1]
Data <- Data[, -1]
Afin de plus facilement repérer les variables et leur modalité, on modifie le jeu de données :
for (j in 1:ncol(Data)) {
if (j < 10) {
Data[, j] <- as.factor(paste("0", j, "-", Data[, j], sep = ""))
} else {
Data[, j] <- as.factor(paste(j, "-", Data[, j], sep = ""))
}
}
Question : Vérifiez la taille du jeu de données et
faites un summary pour contrôler le contenu du jeu de
données Data.
# A FAIRE
Correction :
dim(Data)
[1] 800 10
summary(Data)
odor gill.color stalk.root stalk.color.above.ring
01-n :333 02-b :188 03-?:260 04-w :424
01-f :204 02-p :148 03-b:359 04-p :198
01-y : 69 02-w :113 03-c: 54 04-g : 58
01-s : 63 02-n : 91 03-e:105 04-b : 47
01-a : 42 02-g : 73 03-r: 22 04-n : 39
01-l : 38 02-h : 73 04-o : 17
(Other): 51 (Other):114 (Other): 17
stalk.color.below.ring ring.number ring.type spore.print.color population
05-w :434 06-n: 6 07-e:278 08-w :254 09-a: 30
05-p :194 06-o:734 07-f: 6 08-k :191 09-c: 35
05-g : 48 06-t: 60 07-l:127 08-n :178 09-n: 31
05-b : 47 07-n: 6 08-h :154 09-s:120
05-n : 44 07-p:383 08-r : 9 09-v:413
05-o : 17 08-o : 5 09-y:171
(Other): 16 (Other): 9
habitat
10-d:300
10-g:196
10-l: 93
10-m: 30
10-p:129
10-u: 33
10-w: 19
De même on lit le jeu de données sous Python :
data_full = pd.read_csv('mushroomTP.csv')
target = data_full[['class']]
Datapy = data_full.drop(['class'],axis=1)
Question : Faites une analyse des correspondances multiples avec R. Vous conserverez les coordonnées des projections des individus pour la suite du TP.
library(FactoMineR)
library(factoextra)
resacm <- MCA(....)
fviz_screeplot(resacm, choice = "eigenvalue")
fviz_mca(...)
Correction :
library(FactoMineR)
resacm <- MCA(Data, ncp = 10, graph = FALSE)
fviz_screeplot(resacm, choice = "eigenvalue")
ggplotly(fviz_mca(resacm, geom = "point", geom.var = "text"))
ggplotly(fviz_mca(resacm, geom = "point", geom.var = "text", axes = c(1, 3)))
fviz_mca_ind(resacm, geom = "point", habillage = classesp)
fviz_mca_var(resacm, col.var = "contrib") + scale_color_gradient2(low = "white",
mid = "blue", high = "red", midpoint = 2) + theme_minimal()
fviz_contrib(resacm, choice = "var", axes = 1, top = 10)
Question : Faites également une analyse des correspondances multiples avec Python.
# A COMPLETER
import prince
mca=prince.MCA(....)
mca = mca.fit(Datapy)
Correction :
import prince
mca = prince.MCA(
n_components=10,
n_iter=3,
copy=True,
check_input=True,
engine='auto',
random_state=42
).fit(Datapy)
mca_df = pd.DataFrame({
"Dim1" : mca.row_coordinates(Datapy)[0],
"Dim2" : mca.row_coordinates(Datapy)[1],
"Class": data_full["class"],
"Var1" : mca.column_coordinates(Datapy)[0],
"Var2" : mca.column_coordinates(Datapy)[1]
})
fig=px.scatter(mca_df,x="Dim1",y="Dim2",color="Class")
fig.show()
On se munit également d’une fonction auxiliaire pour visualiser les
résultats du clustering pour la suite. Cette fonction prend en argument
un clustering clust, le jeu de données Data et
un vecteur J contenant les indices des variables d’intérêt.
Il renvoie la fréquence des modalités des variables choisies vis-à-vis
des classes de clust.
# J indice des noms de variables
heatm <- function(clust, Data, J) {
library(dplyr)
Dataaux <- data.frame(id.s = c(1:nrow(Data)), Data)
aux <- cbind(Dataaux, clust)
aux.long <- melt(data.frame(lapply(aux, as.character)), stringsAsFactors = FALSE,
id = c("id.s", "clust"), factorsAsStrings = T)
# Effectifs
aux.long.q <- aux.long %>%
group_by(clust, variable, value) %>%
mutate(count = n_distinct(id.s)) %>%
distinct(clust, variable, value, count)
# avec fréquences
aux.long.p <- aux.long.q %>%
group_by(clust, variable) %>%
mutate(perc = count/sum(count)) %>%
arrange(clust)
Lev <- NULL
for (j in 1:ncol(Data)) Lev <- c(Lev, levels(Data[, j]))
Jaux <- NULL
for (j in 1:length(J)) {
Jaux <- c(Jaux, which(aux.long.p$variable == colnames(Data)[J[j]]))
}
gaux <- ggplot(aux.long.p[Jaux, ], aes(x = clust, y = value)) + geom_tile(aes(fill = perc)) +
scale_fill_gradient2(low = "white", mid = "blue", high = "red") + theme_minimal()
return(list(gaux = gaux, eff = aux.long.q, freq = aux.long.p))
}
Mode <- function(x) {
ux <- unique(x)
ux[which.max(tabulate(match(x, ux)))]
}
barplotClus <- function(clust, Data, J) {
aux.long.p <- heatm(clust, Data, J)$freq
# ux<-unique(aux.long.p$variable)
for (j in J) {
p <- ggplot(aux.long.p[which(aux.long.p$variable == colnames(Data)[j]), ],
aes(x = clust, y = perc, fill = value)) + geom_bar(stat = "identity")
print(p)
}
}
Dans cette partie, nous allons utiliser la méthode des Kmodes introduite par Huang (1998). Rappelons que cette méthode est une extension des Kmeans dans le cas des données qualitatives. Les modifications par rapport aux Kmeans sont
\[ d(\mathbf{x}_i,\mathbf{x}_\ell) = \underset{j=1}{\stackrel{p}{\sum}}\ \mathbb{1}_{x_{ij}\neq x_{\ell j}} \]
\[ \mathbf{m}_k=(m_{k1},\ldots,m_{kp}) \textrm{ avec } m_{kj}= \underset{u_1,\ldots,u_{s_j}}{\mbox{argmax}}\ \underset{i\in\mathcal C_k}{\sum}\ \mathbb{1}_{x_{ij}= u_{s_j}} \]
Question : A l’aide de la fonction
kmodes()de la librairie klaR, déterminez une
classification en \(K=6\) classes des
champignons. Comparez avec la classe de l’espèce. Visualisez la
classification obtenue. Vous pouvez vous aider des fonctions auxiliaires
pour interpréter la classification.
library(klaR)
reskmodes <- kmodes(....)
Correction :
library(klaR)
K <- 6
clkmodes <- kmodes(Data, K, iter.max = 100, weight = FALSE)
table(classesp, clkmodes$cluster)
classesp 1 2 3 4 5 6
e 134 101 29 0 136 0
p 8 38 194 60 32 68
clkmodes$modes
odor gill.color stalk.root stalk.color.above.ring stalk.color.below.ring
1 01-n 02-n 03-b 04-w 05-w
2 01-n 02-p 03-b 04-w 05-w
3 01-y 02-b 03-? 04-w 05-p
4 01-f 02-h 03-b 04-p 05-b
5 01-n 02-p 03-e 04-w 05-w
6 01-f 02-g 03-b 04-b 05-n
ring.number ring.type spore.print.color population habitat
1 06-o 07-p 08-n 09-y 10-d
2 06-o 07-p 08-k 09-v 10-d
3 06-o 07-e 08-w 09-v 10-l
4 06-o 07-l 08-h 09-v 10-d
5 06-o 07-p 08-k 09-s 10-g
6 06-o 07-l 08-h 09-v 10-p
fviz_mca_ind(resacm, habillage = as.factor(clkmodes$cluster), geom = "point") + labs(title = "Kmodes")
fviz_mca_ind(resacm, habillage = as.factor(clkmodes$cluster), geom = "point", axes = c(1,
3)) + labs(title = "Kmodes")
fviz_mca(resacm, habillage = as.factor(clkmodes$cluster), geom.ind = c("point"),
geom.var = c("point", "text"))
# heatm(clust=clkmodes$cluster,Data=Data,J=c(1:10)$gaux
barplotClus(clkmodes$cluster, Data, J = c(1:10))
Question : Pour déterminer le nombre de classes, la méthode du coude peut être utilisée en remplaçant l’inertie intra-classe par le critère “Within Cluster Difference”
\[ WCD(K) = \underset{k=1}{\stackrel{K}{\sum}} \underset{i\in\mathcal C_k}{\sum}\ d(x_i,m_k) \]
où \(d(.,.)\) est l’appariement simple et \(m_k\) est le centre de la classe \(\mathcal C_k\).
Tracez la courbe \(K\mapsto WCD(K)\)
pour déterminer le nombre de classes optimal. Vous pouvez vous aider des
sorties de la fonction kmodes().
# A COMPLETER
WithinDiff <- NULL
Kmax <- 10
Clust <- matrix(0, nrow = nrow(Data), ncol = Kmax)
for (k in 1:Kmax) {
aux <- kmodes(.........)
WithinDiff <- c(WithinDiff, ..........)
Clust[, k] <- aux$cluster
}
auxdf <- data.frame(NbCluster = 1:Kmax, WithinDiff = WithinDiff)
ggplot(auxdf, aes(x = NbCluster, y = WithinDiff)) + geom_point() + geom_line()
Correction :
WithinDiff <- NULL
Kmax <- 10
Clust <- matrix(0, nrow = nrow(Data), ncol = Kmax)
for (k in 1:Kmax) {
aux <- kmodes(Data, k, iter.max = 100, weight = FALSE)
WithinDiff <- c(WithinDiff, sum(aux$withindiff))
Clust[, k] <- aux$cluster
}
auxdf <- data.frame(NbCluster = 1:Kmax, WithinDiff = WithinDiff)
ggplot(auxdf, aes(x = NbCluster, y = WithinDiff)) + geom_point() + geom_line()
Question : Reprenez les questions précédentes avec
Python. Vous pouvez utiliser la fonction KModes.
from kmodes.kmodes import KModes
cost = []
for clustk in range(1, 10):
kmodes = KModes(n_jobs = -1, n_clusters = clustk, init ='Huang',random_state = 0)
kmodes.fit_predict(Datapy);
cost.append(kmodes.cost_);
# A COMPLETER
Correction :
import numpy as np
from kmodes.kmodes import KModes
# Cout
cost = []
for clustk in range(1, 10):
kmodes = KModes(n_jobs = -1, n_clusters = clustk, init ='Huang',random_state = 0)
res=kmodes.fit_predict(Datapy);
cost.append(kmodes.cost_);
# Converting the results into a dataframe and plotting them
df_cost = pd.DataFrame({'Cluster': range(1, 10), 'Cost': cost})
# Data viz
import plotly.express as px
fig = px.line(df_cost, x='Cluster', y='Cost',markers=True)
fig.show()
kmodes = KModes(n_jobs = -1, n_clusters = 6, init ='Huang',random_state = 0);
kmodes.fit(Datapy);
pd.DataFrame(kmodes.cluster_centroids_,columns=Datapy.columns)
odor gill.color stalk.root ... spore.print.color population habitat
0 n n b ... k v d
1 n w ? ... w c g
2 n p e ... n s g
3 f h b ... h v p
4 n p b ... k y d
5 y b ? ... w v l
[6 rows x 10 columns]
pd.crosstab(kmodes.labels_,"freq")
col_0 freq
row_0
0 228
1 62
2 109
3 124
4 80
5 197
Une seconde stratégie est de partir des coordonnées de l’analyse des correspondances multiples et d’utiliser un algorithme plus usuel sur données quantitatives. Dans cette section, on va appliquer l’algorithme des Kmeans.
Question : Appliquez l’algorithme des Kmeans sur les coordonnées de l’ACM. Pour la détermination du nombre de classes, vous pouvez utiliser l’évolution de l’inertie intra-classe et le critère silhouette.
# A COMPLETER
Correction :
set.seed(1234)
Kmax <- 15
kmeansclus <- matrix(0, nrow = nrow(Data), ncol = (Kmax - 1))
Iintra <- NULL
Silhou <- NULL
for (k in 2:Kmax) {
resaux <- kmeans(resacm$ind$coord[, 1:5], centers = k, nstart = 10)
kmeansclus[, (k - 1)] <- resaux$cluster
Iintra <- c(Iintra, resaux$tot.withinss)
aux <- silhouette(resaux$cluster, daisy(resacm$ind$coord))
Silhou <- c(Silhou, mean(aux[, 3]))
}
df <- data.frame(K = 2:Kmax, Iintra = Iintra, Silhou = Silhou)
g1 <- ggplot(df, aes(x = K, y = Iintra)) + geom_line() + geom_point() + xlab("Nombre de classes") +
ylab("Inertie intraclasse")
g2 <- ggplot(df, aes(x = K, y = Silhou)) + geom_line() + geom_point() + xlab("Nombre de classes") +
ylab("Critère Silhouette")
grid.arrange(g1, g2, ncol = 2)
On retient par la suite \(8\) classes.
Question : Etudiez la classification retenue. Comparez avec la classification obtenue précédemment avec les Kmodes et la variable “Class” (champignon comestible ou toxique).
# A COMPLETER
Correction :
Kaux <- 8 - 1
table(kmeansclus[, Kaux])
1 2 3 4 5 6 7 8
17 190 127 45 6 221 177 17
aux <- silhouette(kmeansclus[, Kaux], daisy(resacm$ind$coord))
fviz_silhouette(aux) + theme(plot.title = element_text(size = 9))
cluster size ave.sil.width
1 1 17 0.72
2 2 190 0.76
3 3 127 0.69
4 4 45 0.28
5 5 6 0.92
6 6 221 0.53
7 7 177 0.15
8 8 17 0.57
fviz_mca(resacm, geom.ind = "point", geom.var = c("point", "text"), habillage = as.factor(kmeansclus[,
Kaux]))
fviz_mca_ind(resacm, geom = "point", habillage = as.factor(kmeansclus[, Kaux]), axes = c(1,
2))
fviz_mca_ind(resacm, geom = "point", habillage = as.factor(kmeansclus[, Kaux]), axes = c(1,
3))
# Comparaison avec Kmodes
table(clkmodes$cluster, kmeansclus[, Kaux])
1 2 3 4 5 6 7 8
1 0 0 0 1 0 94 37 10
2 0 0 0 13 3 115 7 1
3 17 190 0 4 2 5 0 5
4 0 0 59 0 1 0 0 0
5 0 0 0 27 0 7 133 1
6 0 0 68 0 0 0 0 0
adjustedRandIndex(clkmodes$cluster, kmeansclus[, Kaux])
[1] 0.5784484
# Comparaison avec comestible/toxique
table(classesp, kmeansclus[, Kaux])
classesp 1 2 3 4 5 6 7 8
e 17 0 0 34 0 186 146 17
p 0 190 127 11 6 35 31 0
adjustedRandIndex(classesp, kmeansclus[, Kaux])
[1] 0.2779427
# heatm(clust=kmeansclus[,9],Data=Data,J=c(1:5))
barplotClus(kmeansclus[, Kaux], Data, J = c(1:10))
Question : Reprenez les questions précédentes avec Python.
# A completer
Correction :
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score
mca = prince.MCA(n_components=5).fit(Datapy)
silhouette_coefficients = []
InertieIntra=[]
for k in range(3, 16):
kmeans = KMeans(n_clusters=k,init = 'random', max_iter = 300, n_init = 10, random_state = 10);
res=kmeans.fit(mca.row_coordinates(Datapy));
score = silhouette_score(mca.row_coordinates(Datapy), kmeans.labels_);
silhouette_coefficients.append(score);
InertieIntra.append(kmeans.inertia_);
aux = pd.DataFrame({"K" : range(3, 16),
"Silhouette" : silhouette_coefficients,
"Intra" : InertieIntra})
fig = px.line(aux, x='K', y='Silhouette',labels={'x':'Number of clusters', 'y':'Silhouette'},markers=True)
fig.show()
fig = px.line(aux, x='K', y='Intra',labels={'x':'Number of clusters', 'y':'Inertie Intra'},markers=True)
fig.show()
from yellowbrick.cluster import SilhouetteVisualizer
kmeansopt = KMeans(6,init = 'k-means++', max_iter = 300, n_init = 10, random_state = 0);
visualizer = SilhouetteVisualizer(kmeansopt, colors='yellowbrick');
visualizer.fit(mca.row_coordinates(Datapy))
SilhouetteVisualizer(ax=<AxesSubplot: >, colors='yellowbrick',
estimator=KMeans(n_clusters=6, random_state=0))
visualizer.show()
kmeanslab=kmeansopt.predict(mca.row_coordinates(Datapy))
pd.crosstab(kmeanslab,"freq")
col_0 freq
row_0
0 127
1 434
2 6
3 19
4 197
5 17
pd.crosstab(kmeanslab,kmodes.labels_)
col_0 0 1 2 3 4 5
row_0
0 0 0 0 124 3 0
1 215 33 109 0 77 0
2 0 6 0 0 0 0
3 0 19 0 0 0 0
4 0 0 0 0 0 197
5 13 4 0 0 0 0
adjusted_rand_score(kmeanslab,kmodes.labels_)
0.515926879468948
pd.crosstab(kmeanslab,data_full["class"])
class e p
row_0
0 0 127
1 359 75
2 0 6
3 19 0
4 5 192
5 17 0
adjusted_rand_score(kmeanslab,data_full["class"])
0.4129657489670915
Question : Quelle classification par CAH pouvez-vous
mettre en place pour traiter des données qualitatives ? Mettez en
application votre proposition avec R et étudiez la classification
retenue. Vous pourrez vous aider des fonctions daisy() de
la librairie cluster, de la fonction hclust(),
…
# A COMPLETER
Correction :
gower.dist <- daisy(Data, metric = c("gower"))
aggl <- hclust(gower.dist, method = "complete")
ggplot(data.frame(K = 1:20, Height = sort(aggl$height, decreasing = T)[1:20]), aes(x = K,
y = Height)) + geom_line() + geom_point()
fviz_dend(aggl, show_labels = FALSE, k = 5)
clustaggl <- cutree(aggl, 5)
fviz_mca(resacm, geom.ind = "point", geom.var = c("point", "text"), habillage = as.factor(clustaggl))
table(clustaggl)
clustaggl
1 2 3 4 5
274 177 216 127 6
table(clustaggl, classesp)
classesp
clustaggl e p
1 199 75
2 177 0
3 24 192
4 0 127
5 0 6
table(clustaggl, clkmodes$cluster)
clustaggl 1 2 3 4 5 6
1 53 54 0 0 167 0
2 89 82 5 0 1 0
3 0 0 216 0 0 0
4 0 0 0 59 0 68
5 0 3 2 1 0 0
barplotClus(clust = clustaggl, Data, J = c(1:10))
Question : Reprenez la question avec python. Vous
pouvez utiliser la librairie gower.
import gower
from scipy.spatial.distance import squareform
from scipy.cluster.hierarchy import dendrogram, linkage
from sklearn.cluster import AgglomerativeClustering
# A COMPLETER
Correction :
import gower
from scipy.spatial.distance import squareform
distance_matrix = gower.gower_matrix(Datapy);
condensed_dist_matrix = squareform(distance_matrix)
from scipy.cluster.hierarchy import dendrogram, linkage
Zd=linkage(condensed_dist_matrix,'complete');
res=dendrogram(Zd, orientation='top',distance_sort='descending',show_leaf_counts=False,no_labels=True);
plt.show()
from sklearn.cluster import AgglomerativeClustering
model = AgglomerativeClustering(n_clusters=9,affinity='precomputed',linkage='complete')
ClustK = model.fit_predict(distance_matrix)
pd.crosstab(ClustK,"freq")
col_0 freq
row_0
0 148
1 196
2 20
3 41
4 162
5 127
6 83
7 6
8 17
pd.crosstab(ClustK,data_full["class"])
class e p
row_0
0 82 66
1 5 191
2 19 1
3 32 9
4 162 0
5 0 127
6 83 0
7 0 6
8 17 0
pd.crosstab(ClustK,kmodes.labels_)
col_0 0 1 2 3 4 5
row_0
0 106 1 41 0 0 0
1 0 0 0 0 0 196
2 0 19 0 0 0 1
3 9 32 0 0 0 0
4 85 0 0 0 77 0
5 0 0 0 124 3 0
6 15 0 68 0 0 0
7 0 6 0 0 0 0
8 13 4 0 0 0 0
adjusted_rand_score(ClustK,kmodes.labels_)
0.6528431966802858
mca_df = pd.DataFrame({
"Dim1" : mca.row_coordinates(Datapy)[0],
"Dim2" : mca.row_coordinates(Datapy)[1],
"Dim3" : mca.row_coordinates(Datapy)[2],
"Clust" : pd.Categorical(ClustK)
})
fig=px.scatter(mca_df,x="Dim1",y="Dim2",color="Clust")
fig.show()
fig=px.scatter(mca_df,x="Dim1",y="Dim3",color="Clust")
fig.show()
Pour cette dernière partie du TP, nous allons utiliser une stratégie de classification via des modèles de mélange. Rappelons que l’on a pour données \(\mathbf{X}=(x_1,\ldots,x_n)\) avec \(x_i\) décrit par \(p\) variables catégorielles, chacune avec \(m_j\) modalités. Dans ce cas de variables qualitatives, on peut considérer des distributions multinomiales par variables et par composantes. Pour écrire la distribution de mélange, on commence par transformer l’écriture des données de la façon suivante :
\[x_i=(x_{i1},\ldots,x_{ip})\rightsquigarrow (x_i^{jh}; j=1,\ldots,p, h=1,\ldots,m_j)\] avec \[ x_i^{\,jh}=\left\{\begin{array}{l l}1 & \textrm{ si } i \textrm{ prend la modalité } h \textrm{ pour la variable }j\\ 0 & \textrm{ sinon.}\end{array}\right. \] Les densités de mélange considérées sont de la forme \[f(.|\theta_K)=\underset{k=1}{\stackrel{K}{\sum}}\pi_k f_k(x_i|\boldsymbol{\alpha}_k)\] avec
Classiquement, on reparamétrise par
\[ \varepsilon_{k}^{\ jh}=\left\{\begin{array}{l l}1 - \alpha_k^{\ jh} & \textrm{ si } a_k^{\ jh}=1\\ \alpha_k^{\ jh} & \textrm{ si } a_k^{\ jh}=0\end{array}\right. \]
Par exemple \((0.3,0.6,0.1) \rightsquigarrow (0,1,0,\ \ 0.3,0.4,0.1)\)
La densité de mélange se réécrit alors \[f(.|\theta_K)=\underset{k=1}{\stackrel{K}{\sum}}\pi_k \underset{j=1}{\stackrel{p}{\prod}}\underset{h=1}{\stackrel{m_j}{\prod}} \left[\left(1-\varepsilon_k^{\,jh}\right)^{a_k^{\ jh}} \left(\varepsilon_k^{\,jh}\right)^{1-a_k^{\ jh}}\right]^{x_i^{jh}}\]
Selon les hypothèses faites sur les \(\varepsilon_{k}^{\ jh}\) et sur les proportions du mélange, on a 10 formes possibles (dans le même esprit que les 24 formes des mélanges gaussiens, cf cours).
Question : A l’aide de la fonction
mixmodCluster de la librairie Rmixmod,
déterminez une classification des données par modèles de mélange. Vous
étudierez en particulier les probabilités conditionnelles
d’appartenance.
library(Rmixmod)
resmixmod <- mixmodCluster(data = Data, nbCluster = 2:15, criterion = c("BIC", "ICL"),
model = mixmodMultinomialModel("Binary_pk_Ekjh"))
# Graphe des critères de sélection
K <- NULL
BIC <- NULL
ICL <- NULL
for (k in 1:length(resmixmod@results)) {
K <- c(K, resmixmod@results[[k]]@nbCluster)
BIC <- c(BIC, resmixmod@results[[k]]@criterionValue[1])
ICL <- c(ICL, resmixmod@results[[k]]@criterionValue[2])
}
## graphique à faire
...
# Etude de la classification retenue
df <- data.frame(proba = apply(resmixmod@bestResult@proba, 1, max), label = as.factor(apply(resmixmod@bestResult@proba,
1, which.max)))
ggplot(df, aes(x = label, y = proba)) + geom_boxplot()
table(resmixmod@bestResult@partition)
## A completer
# Comparaison avec les autres classifications A completer
Correction :
library(Rmixmod)
resmixmod <- mixmodCluster(data = Data, nbCluster = 2:15, criterion = c("BIC", "ICL"),
model = mixmodMultinomialModel("Binary_pk_Ekjh"))
# Graphe des critères de sélection
K <- NULL
BIC <- NULL
ICL <- NULL
for (k in 1:length(resmixmod@results)) {
K <- c(K, resmixmod@results[[k]]@nbCluster)
BIC <- c(BIC, resmixmod@results[[k]]@criterionValue[1])
ICL <- c(ICL, resmixmod@results[[k]]@criterionValue[2])
}
aux <- sort.int(K, index.return = T)
df <- data.frame(K = K[aux$ix], BIC = BIC[aux$ix], ICL = ICL[aux$ix])
df <- melt(df, id = c("K"))
ggplot(df, aes(x = K, y = value, color = variable)) + geom_line()
# Etude de la classification retenue
df <- data.frame(proba = apply(resmixmod@bestResult@proba, 1, max), label = as.factor(apply(resmixmod@bestResult@proba,
1, which.max)))
ggplot(df, aes(x = label, y = proba)) + geom_boxplot()
table(resmixmod@bestResult@partition)
1 2 3 4 5 6 7
69 188 127 126 198 28 64
barplotClus(resmixmod@bestResult@partition, Data, J = c(1:10))
fviz_mca_ind(resacm, habillage = as.factor(resmixmod@bestResult@partition), geom = "point") +
labs(title = "Mixture")
# Comparaison avec Class
table(resmixmod@bestResult@partition, classesp)
classesp
e p
1 69 0
2 0 188
3 0 127
4 77 49
5 181 17
6 24 4
7 49 15
adjustedRandIndex(resmixmod@bestResult@partition, classesp)
[1] 0.2655892
# Comparaison avec Kmeans
table(resmixmod@bestResult@partition, kmeansclus[, Kaux])
1 2 3 4 5 6 7 8
1 0 0 0 0 0 0 69 0
2 0 188 0 0 0 0 0 0
3 0 0 127 0 0 0 0 0
4 0 0 0 0 0 21 105 0
5 0 0 0 0 0 195 3 0
6 17 2 0 4 0 5 0 0
7 0 0 0 41 6 0 0 17
adjustedRandIndex(resmixmod@bestResult@partition, kmeansclus[, Kaux])
[1] 0.8255953